program mainbootstrap
use normalscr2
use matrixfuns
use datahadler
use suplmod
use mardiamodadj
implicit none
integer, parameter:: TT=1000,nn=3,S=10000,sboot=1000
integer,parameter:: nvar=2*nn+5
integer iparam(7),T,n,i,j,i1,i2,si,ind1(TT),ind2(nn),seedvec(S),seedind(2)

double precision:: y(TT,nn),thscale(nvar),rparam(7),a1,a2,c9

double precision:: thlb(nvar),thub(nvar),thtrue(nvar),therr(nvar)&
				 &,th(nvar),llik,grad(nvar),lika,likb,tha(nvar),thb(nvar)

double precision:: mth(S,nvar),suplm(S,4),testmardia(S,3),testjb(S,3)&
                 &,suplmi(sboot,4),testmardiai(sboot,3),testjbi(sboot,3),pvali(sboot)&
				 &,pbootlm(S,4),pbootmardia(S,3),pbootjb(S,3)

double precision:: yi(TT,nn)
double precision, parameter:: pi=3.14159265358979d0
				 
common T,n,y

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
T=TT
n=nn
!!!!!!!!!!True parameters
a1=.1d0
a2=.85d0
c9=.999d0

thtrue(1:n)=vassign(n,.2d0)			!delta
thtrue(n+1:2*n-1)=vassign(n-1,1.0d0)	!c
thtrue(2*n)=0.05d0					!alpha0
thtrue(2*n+1)=0.05d0				!phi0
thtrue(2*n+2)=dasin(dsqrt(a1/c9))    !alpha1
thtrue(2*n+3)=dasin(dsqrt(a2/(c9-a1)))   	!alpha2
thtrue(2*n+4)=dasin(dsqrt(a1/c9))		!phi1
thtrue(2*n+5)=dasin(dsqrt(a2/(c9-a1)))	!phi2
!!!!!!!!! Limits

thlb(1:n)=vassign(n,-1.0d2)			   !delta
thlb(n+1:2*n-1)=vassign(n-1,-1.0d2)		   !c
thlb(2*n)=1.0d-10					   !alpha0
thlb(2*n+1)=1.0d-5					   !phi0
thlb(2*n+2:2*n+5)=vassign(4,0.0d0)

thub(1:n)=vassign(n,1.0d2)			   !delta
thub(n+1:2*n-1)=vassign(n-1,1.0d2)		   !c
thub(2*n)=1.0d2 					   !alpha0
thub(2*n+1)=1.0d2					   !phi0
thub(2*n+2:2*n+5)=vassign(4,2.0d0*pi)

thscale=vassign(nvar,1.0d0)

call readvecint(seedvec,S,"seedsim.txt",1)
call readvecint(seedind,2,"seedind.txt",1)

do si=seedind(1),seedind(2),1  
print*,si
call rnset(seedvec(si))
!!!!!!!!!!Simulation
call gasimul(nvar,T,n,thtrue,y)


!!!!!!!!!!Estimation
call du4inf(iparam,rparam)
iparam(3)=1000
iparam(4)=1000
iparam(5)=1000
!rparam(1)=1.0d-7
!rparam(2)=1.0d-15
call dbcong(gallik,gascr,nvar,thtrue,0,thlb,thub&
          &,thscale,1.0d0,iparam,rparam,th,llik)
mth(si,:)=th

therr=dabs(thub-th)*dabs(th-thlb)
call dsvrgn(nvar,therr,therr)

if (therr(1)>0.0d0) then
	call suplmf(nvar,T,n,y,th,suplm(si,1),suplm(si,2),suplm(si,3),suplm(si,4))
	call mardiajbfadj(nvar,T,n,y,th,testmardia(si,:),testjb(si,:))
	do j=1,sboot,1
		call gasimul(nvar,T,n,th,yi)
	!	call suplmf(nvar,T,n,yi,th,suplmi(j,1),suplmi(j,2),suplmi(j,3),suplmi(j,4))
		call mardiajbfadj2(nvar,T,n,yi,th,testmardiai(j,:),testjbi(j,:))
	end do

!	do i1=1,4,1
!		where (suplmi(:,i1) .gt. suplm(si,i1))
!			pvali=1.0d0
!		elsewhere
!			pvali=0.0d0
!		end where
!			pbootlm(si,i1)=sum(pvali)/dble(sboot)
!	end do

	do i1=1,3,1
		where (testmardiai(:,i1) .gt. testmardia(si,i1))
			pvali=1.0d0
		elsewhere
			pvali=0.0d0
		end where
			pbootmardia(si,i1)=sum(pvali)/dble(sboot)


		where (testjbi(:,i1) .gt. testjb(si,i1))
			pvali=1.0d0
		elsewhere
			pvali=0.0d0
		end where
			pbootjb(si,i1)=sum(pvali)/dble(sboot)
	end do
end if
	call exportvec(vec2(mth(seedind(1):si,:)),"matth.txt",1)
!	call exportrmat(pbootlm(seedind(1):si,:),"pbootlm.txt",1)
	call exportrmat(pbootmardia(seedind(1):si,:),"pbmardia.txt",1)
	call exportrmat(pbootjb(seedind(1):si,:),"pbjb.txt",1)
end do
end program mainbootstrap